//------------------------------------------------------------------------------
// CHOLMOD/Cholesky/t_cholmod_ltsolve_template: template for L'x=b or DL'x=b
//------------------------------------------------------------------------------

// CHOLMOD/Cholesky Module.  Copyright (C) 2005-2023, Timothy A. Davis
// All Rights Reserved.
// SPDX-License-Identifier: LGPL-2.1+

//------------------------------------------------------------------------------

// Template routine to solve L'x=b with unit or non-unit diagonal, or solve
// DL'x=b.
//
// The numeric xtype of L and Y must match.  Y contains b on input and x on
// output, stored in row-form.  Y is nrow-by-n, where nrow must equal 1 for the
// complex or zomplex cases, and nrow <= 4 for the real case.
//
// This file is not compiled separately.  It is included in
// t_cholmod_solve_worker.c instead.  It contains no user-callable routines.
//
// workspace: none
//
// Supports real, complex, and zomplex factors, and any dtype.

// undefine all prior definitions
#undef FORM_NAME
#undef LSOLVE
#undef DIAG

//------------------------------------------------------------------------------
// define the method
//------------------------------------------------------------------------------

#ifdef LL
// LL': solve Lx=b with non-unit diagonal
#define FORM_NAME(prefix,rank) prefix ## ll_ltsolve_ ## rank
#define DIAG

#elif defined (LD)
// LDL': solve LDx=b
#define FORM_NAME(prefix,rank) prefix ## ldl_dltsolve_ ## rank
#define DIAG

#else
// LDL': solve Lx=b with unit diagonal
#define FORM_NAME(prefix,rank) prefix ## ldl_ltsolve_ ## rank

#endif

// LSOLVE(k) defines the name of a routine for an n-by-k right-hand-side.
#define LSOLVE(prefix,rank) FORM_NAME(prefix,rank)

#ifdef REAL

//------------------------------------------------------------------------------
// LSOLVE (1)
//------------------------------------------------------------------------------

// Solve L'x=b, where b has 1 column

static void LSOLVE (PREFIX,1)
(
    cholmod_factor *L,
    Real X [ ]                        // n-by-1 in row form
)
{
    Real *Lx = L->x ;
    Int *Li = L->i ;
    Int *Lp = L->p ;
    Int *Lnz = L->nz ;
    Int j, n = L->n ;

    for (j = n-1 ; j >= 0 ; )
    {
        // get the start, end, and length of column j
        Int p = Lp [j] ;
        Int lnz = Lnz [j] ;
        Int pend = p + lnz ;

        // find a chain of supernodes (up to j, j-1, and j-2)
        if (j < 4 || lnz != Lnz [j-1] - 1 || Li [Lp [j-1]+1] != j)
        {

            //------------------------------------------------------------------
            // solve with a single column of L
            //------------------------------------------------------------------

            Real y = X [j] ;
            #ifdef DIAG
            Real d = Lx [p] ;
            #endif
            #ifdef LD
            y /= d ;
            #endif
            for (p++ ; p < pend ; p++)
            {
                y -= Lx [p] * X [Li [p]] ;
            }
            #ifdef LL
            X [j] = y / d ;
            #else
            X [j] = y ;
            #endif
            j-- ;       // advance to the next column of L

        }
        else if (lnz != Lnz [j-2]-2 || Li [Lp [j-2]+2] != j)
        {

            //------------------------------------------------------------------
            // solve with a supernode of two columns of L
            //------------------------------------------------------------------

            Real y [2], t ;
            Int q = Lp [j-1] ;
            #ifdef DIAG
            Real d [2] ;
            d [0] = Lx [p] ;
            d [1] = Lx [q] ;
            #endif
            t = Lx [q+1] ;
            #ifdef LD
            y [0] = X [j  ] / d [0] ;
            y [1] = X [j-1] / d [1] ;
            #else
            y [0] = X [j  ] ;
            y [1] = X [j-1] ;
            #endif
            for (p++, q += 2 ; p < pend ; p++, q++)
            {
                Int i = Li [p] ;
                y [0] -= Lx [p] * X [i] ;
                y [1] -= Lx [q] * X [i] ;
            }
            #ifdef LL
            y [0] /= d [0] ;
            y [1] = (y [1] - t * y [0]) / d [1] ;
            #else
            y [1] -= t * y [0] ;
            #endif
            X [j  ] = y [0] ;
            X [j-1] = y [1] ;
            j -= 2 ;        // advance to the next column of L

        }
        else
        {

            //------------------------------------------------------------------
            // solve with a supernode of three columns of L
            //------------------------------------------------------------------

            Real y [3], t [3] ;
            Int q = Lp [j-1] ;
            Int r = Lp [j-2] ;
            #ifdef DIAG
            Real d [3] ;
            d [0] = Lx [p] ;
            d [1] = Lx [q] ;
            d [2] = Lx [r] ;
            #endif
            t [0] = Lx [q+1] ;
            t [1] = Lx [r+1] ;
            t [2] = Lx [r+2] ;
            #ifdef LD
            y [0] = X [j]   / d [0] ;
            y [1] = X [j-1] / d [1] ;
            y [2] = X [j-2] / d [2] ;
            #else
            y [0] = X [j] ;
            y [1] = X [j-1] ;
            y [2] = X [j-2] ;
            #endif
            for (p++, q += 2, r += 3 ; p < pend ; p++, q++, r++)
            {
                Int i = Li [p] ;
                y [0] -= Lx [p] * X [i] ;
                y [1] -= Lx [q] * X [i] ;
                y [2] -= Lx [r] * X [i] ;
            }
            #ifdef LL
            y [0] /= d [0] ;
            y [1] = (y [1] - t [0] * y [0]) / d [1] ;
            y [2] = (y [2] - t [2] * y [0] - t [1] * y [1]) / d [2] ;
            #else
            y [1] -= t [0] * y [0] ;
            y [2] -= t [2] * y [0] + t [1] * y [1] ;
            #endif
            X [j-2] = y [2] ;
            X [j-1] = y [1] ;
            X [j  ] = y [0] ;
            j -= 3 ;        // advance to the next column of L
        }
    }
}

//------------------------------------------------------------------------------
// LSOLVE (2)
//------------------------------------------------------------------------------

// Solve L'x=b, where b has 2 columns

static void LSOLVE (PREFIX,2)
(
    cholmod_factor *L,
    Real X [ ][2]         // n-by-2 in row form
)
{
    Real *Lx = L->x ;
    Int *Li = L->i ;
    Int *Lp = L->p ;
    Int *Lnz = L->nz ;
    Int j, n = L->n ;

    for (j = n-1 ; j >= 0 ; )
    {
        // get the start, end, and length of column j
        Int p = Lp [j] ;
        Int lnz = Lnz [j] ;
        Int pend = p + lnz ;

        // find a chain of supernodes (up to j, j-1, and j-2)
        if (j < 4 || lnz != Lnz [j-1] - 1 || Li [Lp [j-1]+1] != j)
        {

            //------------------------------------------------------------------
            // solve with a single column of L
            //------------------------------------------------------------------

            Real y [2] ;
            #ifdef DIAG
            Real d = Lx [p] ;
            #endif
            #ifdef LD
            y [0] = X [j][0] / d ;
            y [1] = X [j][1] / d ;
            #else
            y [0] = X [j][0] ;
            y [1] = X [j][1] ;
            #endif
            for (p++ ; p < pend ; p++)
            {
                Int i = Li [p] ;
                y [0] -= Lx [p] * X [i][0] ;
                y [1] -= Lx [p] * X [i][1] ;
            }
            #ifdef LL
            X [j][0] = y [0] / d ;
            X [j][1] = y [1] / d ;
            #else
            X [j][0] = y [0] ;
            X [j][1] = y [1] ;
            #endif
            j-- ;       // advance to the next column of L

        }
        else if (lnz != Lnz [j-2]-2 || Li [Lp [j-2]+2] != j)
        {

            //------------------------------------------------------------------
            // solve with a supernode of two columns of L
            //------------------------------------------------------------------

            Real y [2][2], t ;
            Int q = Lp [j-1] ;
            #ifdef DIAG
            Real d [2] ;
            d [0] = Lx [p] ;
            d [1] = Lx [q] ;
            #endif
            t = Lx [q+1] ;
            #ifdef LD
            y [0][0] = X [j  ][0] / d [0] ;
            y [0][1] = X [j  ][1] / d [0] ;
            y [1][0] = X [j-1][0] / d [1] ;
            y [1][1] = X [j-1][1] / d [1] ;
            #else
            y [0][0] = X [j  ][0] ;
            y [0][1] = X [j  ][1] ;
            y [1][0] = X [j-1][0] ;
            y [1][1] = X [j-1][1] ;
            #endif
            for (p++, q += 2 ; p < pend ; p++, q++)
            {
                Int i = Li [p] ;
                y [0][0] -= Lx [p] * X [i][0] ;
                y [0][1] -= Lx [p] * X [i][1] ;
                y [1][0] -= Lx [q] * X [i][0] ;
                y [1][1] -= Lx [q] * X [i][1] ;
            }
            #ifdef LL
            y [0][0] /= d [0] ;
            y [0][1] /= d [0] ;
            y [1][0] = (y [1][0] - t * y [0][0]) / d [1] ;
            y [1][1] = (y [1][1] - t * y [0][1]) / d [1] ;
            #else
            y [1][0] -= t * y [0][0] ;
            y [1][1] -= t * y [0][1] ;
            #endif
            X [j  ][0] = y [0][0] ;
            X [j  ][1] = y [0][1] ;
            X [j-1][0] = y [1][0] ;
            X [j-1][1] = y [1][1] ;
            j -= 2 ;        // advance to the next column of L

        }
        else
        {

            //------------------------------------------------------------------
            // solve with a supernode of three columns of L
            //------------------------------------------------------------------

            Real y [3][2], t [3] ;
            Int q = Lp [j-1] ;
            Int r = Lp [j-2] ;
            #ifdef DIAG
            Real d [3] ;
            d [0] = Lx [p] ;
            d [1] = Lx [q] ;
            d [2] = Lx [r] ;
            #endif
            t [0] = Lx [q+1] ;
            t [1] = Lx [r+1] ;
            t [2] = Lx [r+2] ;
            #ifdef LD
            y [0][0] = X [j  ][0] / d [0] ;
            y [0][1] = X [j  ][1] / d [0] ;
            y [1][0] = X [j-1][0] / d [1] ;
            y [1][1] = X [j-1][1] / d [1] ;
            y [2][0] = X [j-2][0] / d [2] ;
            y [2][1] = X [j-2][1] / d [2] ;
            #else
            y [0][0] = X [j  ][0] ;
            y [0][1] = X [j  ][1] ;
            y [1][0] = X [j-1][0] ;
            y [1][1] = X [j-1][1] ;
            y [2][0] = X [j-2][0] ;
            y [2][1] = X [j-2][1] ;
            #endif
            for (p++, q += 2, r += 3 ; p < pend ; p++, q++, r++)
            {
                Int i = Li [p] ;
                y [0][0] -= Lx [p] * X [i][0] ;
                y [0][1] -= Lx [p] * X [i][1] ;
                y [1][0] -= Lx [q] * X [i][0] ;
                y [1][1] -= Lx [q] * X [i][1] ;
                y [2][0] -= Lx [r] * X [i][0] ;
                y [2][1] -= Lx [r] * X [i][1] ;
            }
            #ifdef LL
            y [0][0] /= d [0] ;
            y [0][1] /= d [0] ;
            y [1][0] = (y [1][0] - t [0] * y [0][0]) / d [1] ;
            y [1][1] = (y [1][1] - t [0] * y [0][1]) / d [1] ;
            y [2][0] = (y [2][0] - t [2] * y [0][0] - t [1] * y [1][0]) / d [2];
            y [2][1] = (y [2][1] - t [2] * y [0][1] - t [1] * y [1][1]) / d [2];
            #else
            y [1][0] -= t [0] * y [0][0] ;
            y [1][1] -= t [0] * y [0][1] ;
            y [2][0] -= t [2] * y [0][0] + t [1] * y [1][0] ;
            y [2][1] -= t [2] * y [0][1] + t [1] * y [1][1] ;
            #endif
            X [j  ][0] = y [0][0] ;
            X [j  ][1] = y [0][1] ;
            X [j-1][0] = y [1][0] ;
            X [j-1][1] = y [1][1] ;
            X [j-2][0] = y [2][0] ;
            X [j-2][1] = y [2][1] ;
            j -= 3 ;        // advance to the next column of L
        }
    }
}

//------------------------------------------------------------------------------
// LSOLVE (3)
//------------------------------------------------------------------------------

// Solve L'x=b, where b has 3 columns

static void LSOLVE (PREFIX,3)
(
    cholmod_factor *L,
    Real X [ ][3]                 // n-by-3 in row form
)
{
    Real *Lx = L->x ;
    Int *Li = L->i ;
    Int *Lp = L->p ;
    Int *Lnz = L->nz ;
    Int j, n = L->n ;

    for (j = n-1 ; j >= 0 ; )
    {
        // get the start, end, and length of column j
        Int p = Lp [j] ;
        Int lnz = Lnz [j] ;
        Int pend = p + lnz ;

        // find a chain of supernodes (up to j, j-1, and j-2)
        if (j < 4 || lnz != Lnz [j-1] - 1 || Li [Lp [j-1]+1] != j)
        {

            //------------------------------------------------------------------
            // solve with a single column of L
            //------------------------------------------------------------------

            Real y [3] ;
            #ifdef DIAG
            Real d = Lx [p] ;
            #endif
            #ifdef LD
            y [0] = X [j][0] / d ;
            y [1] = X [j][1] / d ;
            y [2] = X [j][2] / d ;
            #else
            y [0] = X [j][0] ;
            y [1] = X [j][1] ;
            y [2] = X [j][2] ;
            #endif
            for (p++ ; p < pend ; p++)
            {
                Int i = Li [p] ;
                y [0] -= Lx [p] * X [i][0] ;
                y [1] -= Lx [p] * X [i][1] ;
                y [2] -= Lx [p] * X [i][2] ;
            }
            #ifdef LL
            X [j][0] = y [0] / d ;
            X [j][1] = y [1] / d ;
            X [j][2] = y [2] / d ;
            #else
            X [j][0] = y [0] ;
            X [j][1] = y [1] ;
            X [j][2] = y [2] ;
            #endif
            j-- ;       // advance to the next column of L

        }
        else if (lnz != Lnz [j-2]-2 || Li [Lp [j-2]+2] != j)
        {

            //------------------------------------------------------------------
            // solve with a supernode of two columns of L
            //------------------------------------------------------------------

            Real y [2][3], t ;
            Int q = Lp [j-1] ;
            #ifdef DIAG
            Real d [2] ;
            d [0] = Lx [p] ;
            d [1] = Lx [q] ;
            #endif
            t = Lx [q+1] ;
            #ifdef LD
            y [0][0] = X [j  ][0] / d [0] ;
            y [0][1] = X [j  ][1] / d [0] ;
            y [0][2] = X [j  ][2] / d [0] ;
            y [1][0] = X [j-1][0] / d [1] ;
            y [1][1] = X [j-1][1] / d [1] ;
            y [1][2] = X [j-1][2] / d [1] ;
            #else
            y [0][0] = X [j  ][0] ;
            y [0][1] = X [j  ][1] ;
            y [0][2] = X [j  ][2] ;
            y [1][0] = X [j-1][0] ;
            y [1][1] = X [j-1][1] ;
            y [1][2] = X [j-1][2] ;
            #endif
            for (p++, q += 2 ; p < pend ; p++, q++)
            {
                Int i = Li [p] ;
                y [0][0] -= Lx [p] * X [i][0] ;
                y [0][1] -= Lx [p] * X [i][1] ;
                y [0][2] -= Lx [p] * X [i][2] ;
                y [1][0] -= Lx [q] * X [i][0] ;
                y [1][1] -= Lx [q] * X [i][1] ;
                y [1][2] -= Lx [q] * X [i][2] ;
            }
            #ifdef LL
            y [0][0] /= d [0] ;
            y [0][1] /= d [0] ;
            y [0][2] /= d [0] ;
            y [1][0] = (y [1][0] - t * y [0][0]) / d [1] ;
            y [1][1] = (y [1][1] - t * y [0][1]) / d [1] ;
            y [1][2] = (y [1][2] - t * y [0][2]) / d [1] ;
            #else
            y [1][0] -= t * y [0][0] ;
            y [1][1] -= t * y [0][1] ;
            y [1][2] -= t * y [0][2] ;
            #endif
            X [j  ][0] = y [0][0] ;
            X [j  ][1] = y [0][1] ;
            X [j  ][2] = y [0][2] ;
            X [j-1][0] = y [1][0] ;
            X [j-1][1] = y [1][1] ;
            X [j-1][2] = y [1][2] ;
            j -= 2 ;        // advance to the next column of L

        }
        else
        {

            //------------------------------------------------------------------
            // solve with a supernode of three columns of L
            //------------------------------------------------------------------

            Real y [3][3], t [3] ;
            Int q = Lp [j-1] ;
            Int r = Lp [j-2] ;
            #ifdef DIAG
            Real d [3] ;
            d [0] = Lx [p] ;
            d [1] = Lx [q] ;
            d [2] = Lx [r] ;
            #endif
            t [0] = Lx [q+1] ;
            t [1] = Lx [r+1] ;
            t [2] = Lx [r+2] ;
            #ifdef LD
            y [0][0] = X [j  ][0] / d [0] ;
            y [0][1] = X [j  ][1] / d [0] ;
            y [0][2] = X [j  ][2] / d [0] ;
            y [1][0] = X [j-1][0] / d [1] ;
            y [1][1] = X [j-1][1] / d [1] ;
            y [1][2] = X [j-1][2] / d [1] ;
            y [2][0] = X [j-2][0] / d [2] ;
            y [2][1] = X [j-2][1] / d [2] ;
            y [2][2] = X [j-2][2] / d [2] ;
            #else
            y [0][0] = X [j  ][0] ;
            y [0][1] = X [j  ][1] ;
            y [0][2] = X [j  ][2] ;
            y [1][0] = X [j-1][0] ;
            y [1][1] = X [j-1][1] ;
            y [1][2] = X [j-1][2] ;
            y [2][0] = X [j-2][0] ;
            y [2][1] = X [j-2][1] ;
            y [2][2] = X [j-2][2] ;
            #endif
            for (p++, q += 2, r += 3 ; p < pend ; p++, q++, r++)
            {
                Int i = Li [p] ;
                y [0][0] -= Lx [p] * X [i][0] ;
                y [0][1] -= Lx [p] * X [i][1] ;
                y [0][2] -= Lx [p] * X [i][2] ;
                y [1][0] -= Lx [q] * X [i][0] ;
                y [1][1] -= Lx [q] * X [i][1] ;
                y [1][2] -= Lx [q] * X [i][2] ;
                y [2][0] -= Lx [r] * X [i][0] ;
                y [2][1] -= Lx [r] * X [i][1] ;
                y [2][2] -= Lx [r] * X [i][2] ;
            }
            #ifdef LL
            y [0][0] /= d [0] ;
            y [0][1] /= d [0] ;
            y [0][2] /= d [0] ;
            y [1][0] = (y [1][0] - t [0] * y [0][0]) / d [1] ;
            y [1][1] = (y [1][1] - t [0] * y [0][1]) / d [1] ;
            y [1][2] = (y [1][2] - t [0] * y [0][2]) / d [1] ;
            y [2][0] = (y [2][0] - t [2] * y [0][0] - t [1] * y [1][0]) / d [2];
            y [2][1] = (y [2][1] - t [2] * y [0][1] - t [1] * y [1][1]) / d [2];
            y [2][2] = (y [2][2] - t [2] * y [0][2] - t [1] * y [1][2]) / d [2];
            #else
            y [1][0] -= t [0] * y [0][0] ;
            y [1][1] -= t [0] * y [0][1] ;
            y [1][2] -= t [0] * y [0][2] ;
            y [2][0] -= t [2] * y [0][0] + t [1] * y [1][0] ;
            y [2][1] -= t [2] * y [0][1] + t [1] * y [1][1] ;
            y [2][2] -= t [2] * y [0][2] + t [1] * y [1][2] ;
            #endif
            X [j  ][0] = y [0][0] ;
            X [j  ][1] = y [0][1] ;
            X [j  ][2] = y [0][2] ;
            X [j-1][0] = y [1][0] ;
            X [j-1][1] = y [1][1] ;
            X [j-1][2] = y [1][2] ;
            X [j-2][0] = y [2][0] ;
            X [j-2][1] = y [2][1] ;
            X [j-2][2] = y [2][2] ;
            j -= 3 ;        // advance to the next column of L
        }
    }
}

//------------------------------------------------------------------------------
// LSOLVE (4)
//------------------------------------------------------------------------------

// Solve L'x=b, where b has 4 columns

static void LSOLVE (PREFIX,4)
(
    cholmod_factor *L,
    Real X [ ][4]             // n-by-4 in row form
)
{
    Real *Lx = L->x ;
    Int *Li = L->i ;
    Int *Lp = L->p ;
    Int *Lnz = L->nz ;
    Int j, n = L->n ;

    for (j = n-1 ; j >= 0 ; )
    {
        // get the start, end, and length of column j
        Int p = Lp [j] ;
        Int lnz = Lnz [j] ;
        Int pend = p + lnz ;

        // find a chain of supernodes (up to j, j-1, and j-2)
        if (j < 4 || lnz != Lnz [j-1] - 1 || Li [Lp [j-1]+1] != j)
        {

            //------------------------------------------------------------------
            // solve with a single column of L
            //------------------------------------------------------------------

            Real y [4] ;
            #ifdef DIAG
            Real d = Lx [p] ;
            #endif
            #ifdef LD
            y [0] = X [j][0] / d ;
            y [1] = X [j][1] / d ;
            y [2] = X [j][2] / d ;
            y [3] = X [j][3] / d ;
            #else
            y [0] = X [j][0] ;
            y [1] = X [j][1] ;
            y [2] = X [j][2] ;
            y [3] = X [j][3] ;
            #endif
            for (p++ ; p < pend ; p++)
            {
                Int i = Li [p] ;
                y [0] -= Lx [p] * X [i][0] ;
                y [1] -= Lx [p] * X [i][1] ;
                y [2] -= Lx [p] * X [i][2] ;
                y [3] -= Lx [p] * X [i][3] ;
            }
            #ifdef LL
            X [j][0] = y [0] / d ;
            X [j][1] = y [1] / d ;
            X [j][2] = y [2] / d ;
            X [j][3] = y [3] / d ;
            #else
            X [j][0] = y [0] ;
            X [j][1] = y [1] ;
            X [j][2] = y [2] ;
            X [j][3] = y [3] ;
            #endif
            j-- ;       // advance to the next column of L

        }
        else // if (j == 1 || lnz != Lnz [j-2]-2 || Li [Lp [j-2]+2] != j)
        {

            //------------------------------------------------------------------
            // solve with a supernode of two columns of L
            //------------------------------------------------------------------

            Real y [2][4], t ;
            Int q = Lp [j-1] ;
            #ifdef DIAG
            Real d [2] ;
            d [0] = Lx [p] ;
            d [1] = Lx [q] ;
            #endif
            t = Lx [q+1] ;
            #ifdef LD
            y [0][0] = X [j  ][0] / d [0] ;
            y [0][1] = X [j  ][1] / d [0] ;
            y [0][2] = X [j  ][2] / d [0] ;
            y [0][3] = X [j  ][3] / d [0] ;
            y [1][0] = X [j-1][0] / d [1] ;
            y [1][1] = X [j-1][1] / d [1] ;
            y [1][2] = X [j-1][2] / d [1] ;
            y [1][3] = X [j-1][3] / d [1] ;
            #else
            y [0][0] = X [j  ][0] ;
            y [0][1] = X [j  ][1] ;
            y [0][2] = X [j  ][2] ;
            y [0][3] = X [j  ][3] ;
            y [1][0] = X [j-1][0] ;
            y [1][1] = X [j-1][1] ;
            y [1][2] = X [j-1][2] ;
            y [1][3] = X [j-1][3] ;
            #endif
            for (p++, q += 2 ; p < pend ; p++, q++)
            {
                Int i = Li [p] ;
                y [0][0] -= Lx [p] * X [i][0] ;
                y [0][1] -= Lx [p] * X [i][1] ;
                y [0][2] -= Lx [p] * X [i][2] ;
                y [0][3] -= Lx [p] * X [i][3] ;
                y [1][0] -= Lx [q] * X [i][0] ;
                y [1][1] -= Lx [q] * X [i][1] ;
                y [1][2] -= Lx [q] * X [i][2] ;
                y [1][3] -= Lx [q] * X [i][3] ;
            }
            #ifdef LL
            y [0][0] /= d [0] ;
            y [0][1] /= d [0] ;
            y [0][2] /= d [0] ;
            y [0][3] /= d [0] ;
            y [1][0] = (y [1][0] - t * y [0][0]) / d [1] ;
            y [1][1] = (y [1][1] - t * y [0][1]) / d [1] ;
            y [1][2] = (y [1][2] - t * y [0][2]) / d [1] ;
            y [1][3] = (y [1][3] - t * y [0][3]) / d [1] ;
            #else
            y [1][0] -= t * y [0][0] ;
            y [1][1] -= t * y [0][1] ;
            y [1][2] -= t * y [0][2] ;
            y [1][3] -= t * y [0][3] ;
            #endif
            X [j  ][0] = y [0][0] ;
            X [j  ][1] = y [0][1] ;
            X [j  ][2] = y [0][2] ;
            X [j  ][3] = y [0][3] ;
            X [j-1][0] = y [1][0] ;
            X [j-1][1] = y [1][1] ;
            X [j-1][2] = y [1][2] ;
            X [j-1][3] = y [1][3] ;
            j -= 2 ;        // advance to the next column of L
        }

        // NOTE: with 4 right-hand-sides, it suffices to exploit dynamic
        // supernodes of just size 1 and 2.  3-column supernodes are not
        // needed.
    }
}

#endif

//------------------------------------------------------------------------------
// LSOLVE (k)
//------------------------------------------------------------------------------

static void LSOLVE (PREFIX,k)
(
    cholmod_factor *L,
    cholmod_dense *Y,           // nr-by-n where nr is 1 to 4
    cholmod_sparse *Yset        // input pattern Yset
)
{

    #ifdef DIAG
    Real d [1] ;
    #endif
    Real yx [2] ;
    #ifdef ZOMPLEX
    Real yz [1] ;
    Real *Lz = L->z ;
    Real *Xz = Y->z ;
    #endif
    Real *Lx = L->x ;
    Real *Xx = Y->x ;
    Int *Li = L->i ;
    Int *Lp = L->p ;
    Int *Lnz = L->nz ;
    Int n = L->n ;

    ASSERT (L->xtype == Y->xtype) ; // L and Y must have the same xtype
    ASSERT (L->dtype == Y->dtype) ; // L and Y must have the same dtype
    ASSERT (L->n == Y->ncol) ;      // dimensions must match
    ASSERT (Y->nrow == Y->d) ;      // leading dimension of Y = # rows of Y
    ASSERT (L->xtype != CHOLMOD_PATTERN) ;  // L is not symbolic
    ASSERT (!(L->is_super)) ;       // L is simplicial LL' or LDL'

#ifdef REAL
    if (Yset == NULL)
    {

        //----------------------------------------------------------------------
        // real case, no Yset, with 1 to 4 RHS's and dynamic supernodes
        //----------------------------------------------------------------------

        ASSERT (Y->nrow <= 4) ;
        switch (Y->nrow)
        {
            case 1: LSOLVE (PREFIX,1) (L, Y->x) ; break ;
            case 2: LSOLVE (PREFIX,2) (L, Y->x) ; break ;
            case 3: LSOLVE (PREFIX,3) (L, Y->x) ; break ;
            case 4: LSOLVE (PREFIX,4) (L, Y->x) ; break ;
        }

    }
    else
#endif
    {

        //----------------------------------------------------------------------
        // solve a complex linear system or solve with Yset
        //----------------------------------------------------------------------

        ASSERT (Y->nrow == 1) ;

        Int *Ysetp  = (Yset == NULL) ? NULL : Yset->p ;
        Int *Yseti  = (Yset == NULL) ? NULL : Yset->i ;
        Int ysetlen = (Yset == NULL) ? n : Ysetp [1] ;

        for (Int jj = ysetlen-1 ; jj >= 0 ; jj--)
        {

            Int j = Yseti ? Yseti [jj] : jj ;

            // get the start, end, and length of column j
            Int p = Lp [j] ;
            Int lnz = Lnz [j] ;
            Int pend = p + lnz ;

            // y = X [j]
            ASSIGN (yx,yz,0, Xx,Xz,j) ;

            #ifdef DIAG
            // d = Lx [p]
            ASSIGN_REAL (d,0, Lx,p) ;
            #endif
            #ifdef LD
            // y /= d
            DIV_REAL (yx,yz,0, yx,yz,0, d,0) ;
            #endif

            for (p++ ; p < pend ; p++)
            {
                // y -= conj (Lx [p]) * X [Li [p]]
                Int i = Li [p] ;
                MULTSUBCONJ (yx,yz,0, Lx,Lz,p, Xx,Xz,i) ;
            }

            #ifdef LL
            // X [j] = y / d
            DIV_REAL (Xx,Xz,j, yx,yz,0, d,0) ;
            #else
            // X [j] = y
            ASSIGN (Xx,Xz,j, yx,yz,0) ;
            #endif
        }
    }
}

// prepare for the next inclusion of this file in cholmod_solve.c
#undef LL
#undef LD

